## Read some sample ecg data
ecg <- read.table('http://www.indyrad.iupui.edu/public/mmiller3/sample-ecg-1kHz.txt') names(ecg) <- c('t','ecg')

ecg$t <- ecg$t/1000 # convert from ms to s

par(mfrow=c(2,2))

## Plot the ecg:
plot(ecg ~ t, data=ecg, type='l', main='ECG data sampled at 1 kHz', xlab='Time [s]')

## Calculate fft(ecg):
ecg$fft <- fft(ecg$ecg)

## Plot fft(ecg):
#plot(ecg$fft, type='l')

## Plot Mod(fft(ecg)):
plot(Mod(ecg$fft), type='l', log='y', main='FFT of ecg vs index')

## Find the sample period:
delta <- ecg$t[2] - ecg$t[1]

## Calculate the Nyquist frequency:
f.Nyquist <- 1 / 2 / delta

## Calculate the frequencies.  (Since ecg$t is in seconds, delta
## is in seconds, f.Nyquist is in Hz and ecg$freq is in Hz)
## (Note: I may be off by 1 in indexing here ????)

ecg$freq <- f.Nyquist*c(seq(nrow(ecg)/2), -rev(seq(nrow(ecg)/2)))/(nrow(ecg)/2)

## Plot fft vs frequency
plot(Mod(fft) ~ freq, data=ecg, type='l', log='y', main='FFT of ECG vs frequency', xlab='Frequency [Hz]')

## Now let's look at some artificial data: x <- seq(100000)/1000 # pretend we're sampling at 1 kHz

## We'll put in two frequency components, plus a dc offset f1 <- 5 # Hz
f2 <- 2 # Hz
y <- 0.1*sin(2*pi*f1*x) + sin(2*pi*f2*x) + 50 fft.y <- fft(y)
delta <- x[2] - x[1]
f.Nyquist <- 1 / 2 / delta
f <- f.Nyquist*c(seq(length(x)/2), -rev(seq(length(x)/2)))/(length(x)/2)

par(mfrow=c(2,2))
plot(x,y, type='l', xlim=c(0,20))
plot(f, Mod(fft.y), type='l', log='y')

## Now let's zoom in and mark the points were I expect to see peaks: plot(f, Mod(fft.y), type='l', log='y', xlim=c(-10,10)) rug(c(-f1, -f2, 0, f1, f2), col='red', side=3) 